function out=momentsDraws(GG,CC,RR,SDX,ZZ,inSt)
% =========================================================================
% function momentsDraws.m 
% This is a version of momentszselstate designed to be run over draws, for
% which the positions muts be assigned Outside the code. 
%  
% As that function, it computes moments for all the observables in 
% Z, as well as rows of the state vector, differentiating between Moments and IRFs 
% This allows to obtain level IRFS for some variables as well as a larger
% number of states than for the moments. 
% 
%  *Jan 16 2012: Differentiating IRFs from rest of moments Mom*.  
%  *stIn.statePosIRF* can be different from *.statePosMom*. 
%  These are assigned according to the structure *report* which must contain the structures 
%  *report.stateIRF* and *report.stateMom* respectively. 
%  Both should be cells which will be mapped to a row position within .stateNames. 
%  Only report.stateIRF can contain the prefix 'lev'. 
%  If *stIn.statePosIRF* then IRFS will NOT be produced 
%  If *stIn.statePosIRF* then Moments will NOT be produced. 
%  
%% 2. Inputs, inSt structure 
%
%  Here the structure IN is assumed to have all values already assigned for
%  the elements to be recorded over the loop. 
%  See *momentsAssignDefaults.m* for a default assignment function 
% 
%  Fields that are empty will result inSt empty outputs. 
% 
%% 1. Fields
% 
% *.NSimVec*, N to Drop and N to Keep in Simulated Moments
%       [100 | 150 ]  
%
% *.NRep*, Number of times to replicate the simulations 
%       [ 100 ]
%
% *.NCorrel*, Number of Autocorrelations and Autocovariances  
%       [ 4 ]
%
% *.NIRF*, Number of Periods for the Impulse Responses 
%       [ 8 ] 
%
% *.NVarDec*, Number of Decomposition Across Horizons 0,...,NVDec
%       Default is 0 
%
% *.FlagUnitIRF*, IRFs to Unit Impulse or Unit STD
%  
%  IRFs for observables are NOT added 
%
% *.statePosMom* Position of rows in states to extract (rows of GG), for all except the IRFS 
%  
% *.statePosIRF* Position of rows in states to extract (rows of GG), for the
% IRFS Only. Original cell defining this can have the world 'lev' in it,
% which will determine an addition Indicator. 
%
% *.addIRFStates*, Position within extracted states, i.e. rows of the
% states 
% 
%  IRFs for observables are NOT added 
%% 3. Shutting down some computations
% *.NIRF*==0, or *.statePosIRF empty *No IRF
% 
% *IN.NVDEC*==0, No Horizon Decomposition 
% 
% *.IN.NVecSim(2)*=0, No simulated moments 
% 
% *IN.NCorrel*==0, No AC or COV, simulated or Asymptotic 
%
%% 4. OUT structure
%  For Observables (Obs), in which case number of rows, NR, equal to NZ 
%  If inSt.statePosMom non empty, also for Selected States, NR=NS
%
%% I.) Asympotic (Theoretical or Population) Moments 
%
% *.VCV<States or Obs>:*, [NR NR], asymptotic Variance Covariance Matrix for the
% selected states
%
% *.STD<States or Obs>:*, [NR 1], asymptotic STD for the selected states
%
% *.corr<States or Obs>:*, [NR NR NCorr+1], correlations lag 0 through NCorr
%  
% *.varDec<States or Obs>:*, [NR NX], asymptotic variance decomposition. 
%                           
% *.varDecHor<States or Obs>:*, [NR NX NVDEC], Variance decomposition, horizon 1 through 
%                             NVDEC
%
% *Note* State dimension of the IRFS will be different from other variables
% if statePosIRF and statePosMom differ
% 
% *.IRF<States or Obs>:*, [NR NX NIRF+1], Impulse Responses 
%
%% Simulated (Sample) Moments 
% 
% .STDSim<States or Obs>      [NR NPAGE]
%
% .correlSim<States or Obs>   [NR NR NPAGE+1] 
% 
% Alejandro Justiniano (C)
%
% Created January 2nd 2012
%

%% 1. Check dimension and default values 

if nargin < 6; error('Need 6 inputs');end 
NY=size(GG,1);
NX=size(SDX,1);
NZ=size(ZZ,1);
if isfield(inSt,'statePosMom') &&~isempty(inSt.statePosMom)
    NS=length(inSt.statePosMom);
else
    NS =0;
end
    

%% 2. Theoretical (Monthly) Population Moments

%% 2.1 VCV [NR NR] Theoretical Variances 
pshat=lyapunov_symm(GG,RR*(SDX'*SDX)*RR');
if NS > 0
    out.VCVStates=pshat(inSt.statePosMom,inSt.statePosMom);
else
    out.VCVStates=[];
end
out.VCVObs=ZZ*pshat*ZZ';

%% 2.2. STD [NR  1] Theoretical Standard Deviations 
out.STDObs=sqrt(diag(out.VCVObs));
if NS > 0
    out.STDStates=sqrt(diag(out.VCVStates));
    denos=out.STDStates*out.STDStates'; 
else
    out.STDStates=[];
end
deno=out.STDObs*out.STDObs';

%% 2.3. correl [NR NR NPER+1] Theoretical Correlations 
out.corrObs= zeros(NZ,NZ,inSt.NCorrel+1);
out.covObs   = zeros(NZ,NZ,inSt.NCorrel+1);
if NS > 9
    out.corrStates= zeros(NS,NS,inSt.NCorrel+1);
    out.covObs= zeros(NZ,NZ,inSt.NCorrel+1);
else
    out.corrStates=[];
    out.covObs=[]; 
end
%% 2.4. varDec [NR NX] Stationary Variance decomposition
out.varDecObs=zeros(NZ,NX);
if NS > 0
    out.varDecStates=zeros(NS,NX);
else
    out.varDecStates=[];
end
pall=zeros(NY,NY,NX); 

%% 3. Computation Stationary Variance Decomposition 
for ii=1:NX;    
    pall(:,:,ii)=real(disclyap_fast(GG,RR*((SDX(ii,:)' )*SDX(ii,:))*RR'));
    out.varDecObs(:,ii)=diag(squeeze(ZZ*pall(:,:,ii)*ZZ'))./(diag(out.VCVObs));    
    if NS > 0
        out.varDecStates(:,ii)=diag( pall(inSt.statePosMom,inSt.statePosMom,ii) )./(diag(out.VCVStates));
    end            
end

%% 4. Computation of Covariances and Autocorrelations 
out.corrObs(:,:,1)=squeeze(out.VCVObs(:,:,1))./deno;
out.covObs(:,:,1)=out.VCVObs(:,:,1); 
if NS > 0
    out.covStates(:,:,1)=out.VCVStates(:,:,1); 
    out.corrStates(:,:,1)=squeeze(out.VCVStates(:,:,1))./denos;
    
end
if inSt.NCorrel > 0
    Gprod=eye(NY);
    for ii=1:inSt.NCorrel;
        Gprod=Gprod*GG;
        temp=Gprod*pshat;
        if NS > 0
            out.covStates(:,:,ii+1)=temp(inSt.statePosMom,inSt.statePosMom);
            out.corrStates(:,:,ii+1 )=squeeze(out.covStates(:,:,ii+1))./denos;
        end
        out.covObs(:,:,ii+1)=ZZ*temp*ZZ';
        out.corrObs(:,:,ii+1 )=squeeze(out.covObs(:,:,ii+1))./deno;
    end
end 
%% 5. Computation of Impulse Responses     
if inSt.flagUnitIRF~=1 && inSt.flagUnitIRF~=0
    inSt.flagUnitIRF=0;
end
if inSt.NIRF~=0
    out.IRFObs=zeros(NZ,NX,inSt.NIRF+1);
    if isfield(inSt,'statePosIRF') && ~isempty(inSt.statePosIRF);
        NSIrf=length( inSt.statePosIRF );
        out.IRFStates=zeros(NSIrf,NX,inSt.NIRF+1);
    else
        out.IRFStates=[];
        NSIrf=0; 
    end
    for ii=1:NX
        if inSt.flagUnitIRF==0
            resp=SDX(ii,ii)*RR(:,ii);
        else
            resp=RR(:,ii);
        end
        if NSIrf> 0
            out.IRFStates(:,ii,1)=resp(inSt.statePosIRF);
        end
        out.IRFObs(:,ii,1)=ZZ*resp;
        for mm=2:(inSt.NIRF+1);
            resp=GG*resp;
            if NS > 0
                out.IRFStates(:,ii,mm)=resp(inSt.statePosIRF);
            end
            out.IRFObs(:,ii,mm)=ZZ*resp;
        end
        
        % IRFS for observables NO Longer added 
        %if ~isempty(inSt.addIRFObs);
        %    out.IRFObs(inSt.addIRFObs,ii,:)=cumsum(out.IRFObs(inSt.addIRFObs,ii,:),3);
        %end
        if NSIrf > 0
            if ~isempty(inSt.addStateIRF);
                out.IRFStates(inSt.addStateIRF,ii,:)=cumsum(out.IRFStates(inSt.addStateIRF,ii,:),3);
            end
        end
        
    end
else
    out.IRFObs=[];
    out.IRFStates=[];
end

%% 6. Variance decompositions at horizon 0 through NVDEC
if inSt.NVarDec==0;
    out.varDecHorObs=[];
    out.varDecHorStates=[]; 
else
    tempOut=varDecHor(GG,RR,ZZ,SDX,inSt.NVarDec);
    out.varDecHorObs=tempOut.varDecHorObs; 
    if NS > 0
        out.varDecHorStates=tempOut.varDecHorStates(inSt.statePosMom,:,:);
    end
end

%% II. Simulated (Monthly) Sample Moments

%% II.1 Storage
if inSt.NSimVec(2)==0 
     out.STDSimObs=[]; 
     out.STDSimStates=[]; 
     out.corrSimObs=[];
     out.corrSimStates=[]; 
     return
end
nobs=inSt.NSimVec(1)+inSt.NSimVec(2);
out.STDSimObs=zeros(NZ,inSt.NRep);
if NS > 0
    out.STDSimStates=zeros(NS,inSt.NRep);
end
if inSt.NCorrel > 0
    if NS > 0
        out.corrSimStates=zeros(NS,NS,inSt.NCorrel+1,inSt.NRep);
    else
        out.corrSimStates=[];
    end
    out.corrSimObs=zeros(NZ,NZ,inSt.NCorrel+1,inSt.NRep);
else 
    out.corrSimStates=[]; 
    out.corrSimObs   =[]; 
end 

%% II.2 Simulation 
for hh=1:inSt.NRep;
    shocks=(randn(nobs,NX)*SDX)';
    ssim=zeros(size(GG,1),nobs);
    ysim=zeros(size(ZZ,1),nobs);
    ssim(:,1)=RR*shocks(:,1);
    ysim(:,1)=ZZ*(ssim(:,1) + CC);
    for ii=2:nobs;
        ssim(:,ii)=GG*ssim(:,ii-1)+RR*shocks(:,ii);
        ysim(:,ii)=ZZ*(ssim(:,ii) + CC);
    end
    ysim=ysim(:,inSt.NSimVec(1)+1:end)';
    ssim=ssim(:,inSt.NSimVec(1)+1:end)';
    out.STDSimObs(:,hh)=(std(ysim)');
    if NS > 0
        out.STDSimStates(:,hh)=(std(ssim(:,inSt.statePosMom))');
    end    
    if inSt.NCorrel > 0
        [out.corrSimObs(:,:,:,hh)]=accov(ysim,(0:inSt.NCorrel));
        if NS > 0
            [out.corrSimStates(:,:,:,hh)]=accov(ssim(:,inSt.statePosMom),(0:inSt.NCorrel));
        end
    end
end

end
